Energy spectrum of harmonically trapped two-component Fermi gases: Three- and 

Four-Particle Problem 
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Trapped two-component Fermi gases allow for the investigation of the so-called BCS-BEC 
crossover by tuning the interspecies atom-atom s-wave scattering length scattering a' aa ' from at- 
tractive to repulsive, including vanishing and infinitely large values. Here, we numerically determine 
the energy spectrum of the equal-mass spin-balanced four-fermion system — the smallest few-particle 
system that exhibits BCS-BEC crossover-like behavior — as a function of a' aa ' using the stochastic 
variational approach. For comparative purposes, we also treat the two- and three-particle systems. 
States with vanishing and finite total angular momentum as well as with natural and unnatural parity 
are considered. In addition, the energy spectrum of weakly-attractive and weakly-repulsive gases is 
characterized by employing a perturbative framework that utilizes hyperspherical coordinates. The 
hyperspherical coordinate approach allows for the straightforward assignment of quantum numbers 
and furthermore provides great insights into the strongly-interacting unitary regime. 

PACS numbers: 03.75.Ss,05.30.Fk,34.50.-s 



I. INTRODUCTION 



Few-body systems show rich behaviors that range from 
the realization of highly-correlated states to weakly- 
bound Borromean states, and have long been of great 
interest to chemists as well as nuclear and atomic physi- 
cists. To date, the determination of the entire energy 
spectrum, or parts thereof, of small bosonic or fermionic 
systems consisting of four or more constituents remains 
a challenge despite the ever increasing computational re- 
sources. Recently, significant progress has been made in 
the theoretical characterization of weakly-bound bosonic 
tretramers [H41]- In particular, for each Efimov trimer 
there exist two tetramer states 0, Q, which dissoci- 
ate into four free bosons at critical negative scattering 
lengths. The ratio of the scattering length at which the 
trimer state becomes unbound and that at which the first 
or second tetramer state becomes unbound has been pre- 
dicted to be universal Recently, this prediction has 
been confirmed by loss rate measurements on the nega- 
tive scattering length side Q. Working with an atomic 
Cs sample at temperatures just above the transition tem- 
perature to quantum degeneracy, the Innsbruck group Q 
was able to observe enhanced losses at magnetic field 
strengths that correspond quite well to the theoretically 
predicted scattering length ratios [f|. By now, several 
groups have reported experimental evidence for univer- 
sal four-boson physics 



While the universal properties of few-boson systems in- 
teracting through short-range potentials depend on two 
atomic physics parameters, i.e., the two-body s-wave 
scattering length cS aa ^ and a three-body parameter (see, 
e.g., Ref. the universal properties of dilute equal- 

mass two-component Fermi gases interacting through 
short-range potentials with interspecies s-wave interac- 
tions depend only on the s-wave scattering length fill — 
[28j . Experimentally, small two-component Fermi gases 



can be realized by loading a deep three-dimensional op- 
tical lattice with a deterministic number of atoms per 
lattice site [29rl3lj . If the tunneling between neighboring 
sites is negligible, each lattice site can be treated as an 
independent, approximately harmonically confined few- 
fermion system. 

This paper determines and characterizes the energy 
spectrum of three- and four-particle equal-mass two- 
component Fermi gases as a function of the s-wave scat- 
tering length under spherically symmetric harmonic con- 
finement. For this confining geometry, the total angular 
momentum L to t and the total parity II tot are good quan- 
tum numbers throughout the entire BCS-BEC crossover. 
The three-fermion spectra have been discussed previ- 
ously [H-dll and are included here primarily for illus- 
trative and comparative purposes. Our calculations fol- 
low two distinctly different avenues. On the one hand, 
we numerically determine the energy spectrum of few- 
fcrmion systems throughout the entire crossover. For the 
four-fermion system, we employ the stochastic variational 
app roach (35H4l| . In contrast to previous studies [35|, [40l — 
|43|, we utilize basis functions with well-defined angular 
momentum and parity and determine the eigenenergies 
for a range of angular momenta. On the other hand, 
we determine the eigenspectrum semi-analytically within 
first order degenerate perturbation theory. While neces- 
sarily limited to small |a^°°^|, this approach allows for the 
classification of a large portion of the energy spectrum in 
terms of appropriate quantum numbers. To characterize 
the energy spectrum in the weakly- interacting regime, we 
employ hyperspherical coordinates [20L |43 - |51| and write 
the non-interacting wave functions in the relative coor- 
dinates as a product of a hyperangular channel function 
and a hyperradial weight function. The eigenenergies of 
the non-interacting system have, in general, large degen- 
eracies, which are partially lifted by the two-body inter- 
actions. The energy splittings can, to leading order, be 
calculated perturbatively. Compared to calculations that 
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utilize Cartesian single particle coordinates, one distinct 
advantage of the hyperspherical approach is that cer- 
tain features carry over, with some modifications, to the 
strongly-interacting unitary regime [20|, IH, HH, Hq|> HH ■ 
Our numerically determined spectra at unitarity can thus 
be interpreted within the hyperspherical framework. 

The remainder of this manuscript is organized as 
follows. Section III Al introduces the system Hamilto- 
nian under study and provides other background in- 
formation. Section III Bl discusses the hyperspherical 
framework and its implications for the non-interacting, 
weakly-interacting and strongly-interacting three- and 
four-fermion systems. The numerical basis set type ex- 
pansion approaches for the three- and four-fermion prob- 
lems are discussed in Sec. Ill Cl and Sec. Ill Dl respectively. 
Section Hni summarizes our numerical and semi-analytical 
perturbative results for the three- and four-fermion sys- 
tems. We discuss the degeneracies and quantum numbers 
of the energy levels throughout the BCS-BEC crossover. 
Furthermore, we characterize the energy spectrum at uni- 
tarity and present a simple model that predicts the en- 
ergy spectrum of the three-fermion system and a subset 
of the energy spectrum of the four-fermion system at uni- 
tarity. Lastly, Sec. II VI summarizes our main results. 



II. THEORETICAL BACKGROUND 

This section introduces the system Hamiltonian and 
discusses our approaches to determining the eigenspec- 
trum of equal-mass two-component Fermi gases pertur- 
batively and numerically. 



The potential V5 nt accounts for the short-range two-body 
interactions Vtb between unlike atoms, 



N-, 



N 



E 

j=l k=N t +l 



Vtb(rj-f k ), 



(4) 



where the number of spin-up atoms and the number 
Ni of spin-down atoms add up to the total number of 
atoms, i.e., Nf + N± = N. For spin-imbalanced systems, 
TV-f- denotes the number of atoms of the majority species 
and Ni that of the minority species. Throughout, we as- 
sume that the two-body potential Vtb is characterized by 
the s-wave atom-atom scattering length a^ aa ^ and possi- 
bly a range parameter ro . The different functional forms 
of Vtb employed in our calculations are discussed below. 
The goal of this paper is to determine and interpret the 
eigenenergies E(Nf,Ni) of the Hamiltonian H, Eq. (p}. 

If the atom-atom scattering length a^ aa ^ is negative 
and small in absolute value, i.e., |a^ aa -*| -C aiio, where aho 
denotes the oscillator length associated with the atom 
mass to, 



h 



(5) 



then the Fermi system behaves like a weakly-attractive 
atomic gas. In this case, the energy shifts due to the in- 
teractions can be described, to leading order, within first 
order degenerate perturbation theory that treats H m as 
the unperturbed Hamiltonian and VS n t as the perturba- 
tion (35l l4lj . It is then convenient to parametrize the 
two-bod y p otential Vtbifjk) by Fermi's pseudo-potential 
VFifjk) [HI, 



A. System Hamiltonian and other background 
information 



We consider small equal-mass two-component Fermi 
gases under external harmonic confinement consisting of 
N atoms with mass to and position vectors fj , measured 
with respect to the center of the trap. Our model Hamil- 
tonian H reads 

H = H ai + V int (f 1 ,--- ,r N ), (1) 

where the non-interacting Hamiltonian H nl is given by 

N 



3=1 



(2) 



and the external spherically symmetric harmonic confin- 
ing potential Vtrap is characterized by the angular trap- 
ping frequency u, 



Vtrap (r 3 ) = r j 



(3) 



vH^fe) 



Anh 2 



laa) 



(6) 



which allows for an analytical evaluation of the matrix 
elements and, if employed within first order perturbation 
theory, does not lead to divergencies. In general, multiple 
eigenfunctions ipj 1 of the non-interacting atomic system 
are degenerate so that the first order energy shifts £W 
are obtained by solving the determinantal equation 



det (y int - = 0, 



(7) 



where / denotes the identity matrix and the matrix ele- 
ments (V_ iat )jk are given by 



(V int ) jk = (VflVintIC 



(8) 



Here, j and k run from 1 to g n i, where <7 n i denotes the 
degeneracy of the eigenenergy E n i(Nf, N±) under consid- 
eration. 

The possibly most direct approach for constructing the 
eigenfunctions tfjj 1 is to write the tpj 1 as a product of 
two determinants, one for the spin-up atoms and one 
for the spin-down atoms. The determinants themselves 
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are constructed from the single-particle wave functions 
€ P p ,l p , mp (f P ), P = !,-••, iV t or V = N t + 1, • • • , A, which 
arc cigenfunctions of the single-particle harmonic oscilla- 
tor Hamiltonian H p p , 



H. 



SP 



h 2 „ 

2m rp 



Vtrap^p). 



(9) 



Above, n p , Z p and m p denote the single particle radial, 
orbital angular momentum and projection quantum num- 
bers, respectively. For a given energy E^Nf, AjJ of the 
non-interacting many-body system, the single-particle 
wave functions l m (f p ) have to be chosen such that 
their eigenenergies obey the constraint 

E ni (N t ,N l ) = J2(^P + 



JV 



53 ( 2n p + 1 p + ? ) fiw, ( 10 ) 

p=iV t +l 



with the additional restriction that the sets of quantum 
numbers (n p ,lp,m p ) and (n q ,l q ,m q ) differ by at least 
one entry for p ^ q, where p, q — 1, ■ ■ • , N-f or p, q = 
Nf + 1,- • • , A. Following this approach, the first order 
energy shifts EW of the energetically lowest lying gas- 
like state for systems with Nf — N± = 0, ±1 and up to 
A = 20 atoms have been calculated (35| . 

Although in principle straightforward, the outlined 
construction of the wave functions of the non-interacting 
Fermi gas and their use in evaluating the energy shifts 
E^ has several disadvantages. The number of degener- 
ate states of the non-interacting system increases rapidly 
with increasing energy. For example, for the three- 
particle system with (A-f,Aj,) — (2,1), the lowest four 
energies £ , ni (2, 1) = llfej/2, 13^/2, 15^/2 and 17fiw/2 
of the non-interacting system have degeneracies g n i = 3, 
18, 73 and 228, and the determination of the energy shifts 
thus requires the construction and diagonalization of in- 
creasingly large interaction potential matrices V_ int . Fur- 
thermore, the outlined construction includes center-of- 
mass excitations and does not take advantage of the fact 
that the total angular momentum L to t , the corresponding 
z-projection and the parity II tot of the system are good 
quantum numbers. Lastly, the anti-symmetrization is ac- 
complished through the use of determinants, leading to 
Nf\ x A^! terms for each iffi. 

This paper pursues an alternative approach and writes 
the non-interacting wave functions in terms of hyper- 
spherical coordinates 0, l44U5lj ]. This approach sepa- 
rates off the center-of-mass degrees of freedom, treats 
one angular momentum at a time and ensures the proper 
anti-symmetry of the wavefunction by utilizing angular 
momentum algebra. Using the wave functions of the non- 
interacting atomic Fermi gas, written in terms of hyper- 
spherical coordinates, we are able to determine the first- 
order energy shifts for a large portion of the spectrum 
of weakly-interacting equal-mass two-component atomic 



Fermi gases with A = 3 and 4 semi-analytically (see 
Sec. US] and Sec. Hill). 

When a,( aa ) is positive and small (a^ aa ^ <C <2ho) ; di- 
atomic bosonic molecules can form and, if this happens, 
the Fermi system behaves like a weakly-repulsive molec- 
ular Bose gas. In this limit, the dimers or diatomic 
molecules can to a good approximation be treated as 
bosonic point particles with mass 2m pH , l35l . l4ll l54l |55| 
and internal energy -Edimci-; as detailed below, this in- 
ternal energy accounts for the presence of the external 
confinement. The effective model systems for A = 3 and 
A = 4 then consist of two particles, an atom and a dimcr 
in the three-particle case and two dimers in the four- 
particle case [35|, [54Tj56| . Separating off the center-of- 
mass motion, the dynamics are governed by the relative 
effective Hamiltonian H eS , 



h 2 1 
2fiW r 2^ 



(11) 



where k stands for ad (atom-dimer) and dd (dimer-dimer) 
for the three- and four-fermion systems, respectively, and 
the position vector r denotes the atom-dimer and dimer- 
dimer distance vector for the three- and four-fermion sys- 
tems, respectively. The reduced masses // ad ) and ^S dd ^ of 
the atom-dimer and dimer-dimer systems are defined as 
/j( ad ) = 2m/3 and p^ dd ^ = m. In this model, the atom- 
molecule and molecule-molecule interactions are conve- 
niently described through Fermi's regularized pseudo- 
potential V^t f) H3, 



-a w S(r)—r, 



dr 



(12) 



with effective atom-dimer and dimer-dimer scattering 
lengths and a^ dd \ respectively [H [U |H , 



and M, M, m 



j(od) ~ i.i790662349a (aa) 



a {dd) « 0.608a( aa) . 



(13) 



(14) 



The s-wave (l e s = 0) eigenenergies E e g of H eS are readily 
obtained by solving the transcendental equation 61] 



q( fc ) 
,(*) 



r (-Ems, 

L V 2hio 



or I -Ecff 

Z1 V 2hui 



(15) 



where 



(16) 



Since Fermi's regularized zero-range potential Vp k ^ c<r only 

acts at r — 0, the eigenstates of H cS with non- vanishing 
angular momentum l e g do not feel the interaction and 
the corresponding eigenenergies coincide with those of 
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relative angular momentum L IC \ of the three- and four- 
fcrmion systems is and the states have natural par- 
ity, i.e., IT rc i = +1. If E c ff is taken to be an eigenen- 
ergy of H with finite angular momentum Z e g, i.e., 
E c s — (2n e g + I c b + 2>/2)hw with n e g — 0, 1, • • • , then 
the total relative angular momentum of the three- and 
four-fermion systems is L re i = i e g and the states have, as 
above, natural parity, i.e., n rc i = ( — l) Llcl . These obser- 
vations are used below to interpret Figs. [51 and [TOl 



(aa) 



a, / a 

ho 



FIG. 1: Relative s-wave energies E TC i(l, 1) of the trapped 
atom-atom system, obtained by solving Eq. (|15[) for k = aa, 
as a function of aho/a' aa '. The lowest s-wave eigenenergy 
(solid line) is referred to as -Edimor in the text. Here, Eho 
denotes the harmonic oscillator energy, Eh = Tac- 



tile non-interacting system. Figure Q] shows the relative 
eigenenergies B re i(l,l) obtained by solving Eq. (fl~5l) for 
k = aa as a function of a\ lo /a^ aa ^ 1 [in this case, E e g = 
E lel (l,l) and = m/2]. From Eq. JI5J) one finds 
at unitarity B U nit,rel(l> 1) = (2n e g + l/2)Sw for s-wave 
states. The eigenenergies E e g of the effective atom-dimcr 
and dimcr-dimer systems can be obtained from Fig. [T] by 
scaling the horizontal axis appropriately. 

Within the effective two-particle model, the relative 
energies -E re i(2,l) and -E re i(2,2) of the three- and four- 
fermion systems are given by E YC \(2, 1) = E c g + -Edimcr 
and -E re i(2,2) = E cS + 2£ dimor , where the second term 
on the right hand sides accounts for the internal molec- 
ular binding energy of the dimer(s) in the presence of 
the trap [HI, IH| . -Edimer is given by the lowest energy 
solution of Eq. (fT5)) with k = aa and ^- aa ^ = m/2 (see 
solid line in Fig. [T]). For a^ aa ^ <C aho, the size of the 
dimer — given to first order by a( aa ) — is much smaller 
than the trap size and -Edimer approaches the free-space 
result £&- r 



-h 2 /[m(a ( - aa '>) 2 ]. As a (aa) increases, the 
role of the confinement becomes increasingly more impor- 
tant and the lowest energy solution of Eq. (fT5|) starts to 
deviate from E^^ eI . Expanding Eq. (fl~5j) about the non- 
interacting energies (2n c g + 3/2)Huj, the s-wave energies 
E e g can be approximated as [6l| 



2n r 



3/2 



2r(n eff + 3/2) 



v^rKff + i)r(3/2) gW 

ho, fi 




hu (17) 



for small la^l- Alternatively, this result can be obtained 
by treating the atom-dimer and dimer-dimer interactions 
within first order perturbation theory. 

Lastly, we discuss the angular momentum L IC \ of the 
three- and four-fermion systems implied by the effective 
model Hamiltonian H cS . If E c g is taken to be one of 
the positive energy solutions of Eq. (jT5j) , then the total 



B. Hyperspherical Coordinate Approach 

The hyperspherical framework serves two distinct pur- 
poses in this paper. It allows for (i) the construction 
of non-interacting wave functions with good quantum 
numbers and (ii) the classification of the energy spec- 
trum at unitarity. This section first treats the non- 
interacting Fermi gas using hyperspherical coordinates 
and then reviews how the formalism, with some modifi- 
cations, carries over to the infinitely strongly-interacting 
unitary Fermi gas. 

To construct the eigenfunctions of the non-interacting 
Fermi gas, we write the many-body Hamiltonian H ni , 
Eq. 0, in hyperspherical coordinates [2(1 llil - tol"} . We 
first separate off the center-of-mass vector R cm and then 
divide the remaining 3N — 3 coordinates into 3./V — 4 hy- 
perangles, collectively denoted by f2 (see below for their 
definition), and the hyperradius R, 



R 2 



1 - 

-5>-i? cm ) 2 



(18) 



Using these coordinates, the Hamiltonian H m can be 
written as 



H m = H cln — 



h 2 f d 2 

2M \dR 2 



3JV-4 a 



R dR 



A 2 



1 



2MR 2 2 

where the center-of-mass Hamiltonian H cm 



Mu 2 R 2 , (19) 



— V 2 - 

2M «c 



1 



Mto 2 R 2 CTa 



is given by 
(20) 



and M denotes the total mass of the system, i.e., M = 
Nm. In Eq. (|T9")) , A denotes the so-called grand angular 
momentum operator [44[ that accounts for the kinetic 
energy associated with the hyperangles il. The eigen- 
functions of the Hamiltonian H nl separate (see, e.g., 
Refs. Hi H, 



G 



N cm ,L om ,M c 



,r N ) = 

XRcm)F q , x (R)$x, x (ti)- 



(21) 



Here, the center-of-mass functions GN cm ,L ota ,M cm (Rem) 
are eigenfunctions of H cm , i.e., three-dimensional har- 
monic oscillator functions in the center-of-mass vector 
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R cm , with eigenenergies E cm = (2N cm + A cm + 3/2)hu, 
where N cm = 0, 1, • • • , A cm = 0,1,- •• and M cm 

/.. — —Lcm + 1, • • ■ , A cm . The hyperspherical harmon- 
ics $a.x(^)i or so-called channel functions, are eigenfunc- 
tions of the operator A 2 [4J|, 



A 2 $ AjX (r2) = h 2 X(X + SN- 5)$ A , x (fi), 



(22) 



where A can take the values 0, 1, 2, Th e q uantum 
number \ denotes the degeneracy for each A [44J, 



X 



(3N + 2A- 5)(3iV + A -6)! 

A!(3Af-5) ' 



(23) 



In deriving Eq. (|23[) . no symmetry constraints have 
been enforced. Below, we discuss the construction 
of the hyperspherical harmonics and the reduction of 
the degeneracy x due to symmetry constraints for the 
three- and four-fermion systems. Since the center-of- 
mass coordinates and the hyperradius are unchanged un- 
der the exchange of fj and ffe (j, k = 1, ■•• ,JV), the 
symmetry-constraints only effect the $x )X (fi!) an d neither 

GjV cm ,L cm ,M cm (A cm ) nOr F g< x(R). 

Plugging Eq. (|2"Tj) into the Schrodinger equation 
H ni tfj ni = A^A^, A i )?/> ni and dividing out the center- 
of-mass and hyperangular contributions, we obtain an 
effective hyperradial Schrodinger equation [2(], H3| , 



d 2 , ft 2 A' m (K ni + 1) , 1 



2M <9A 2 



2Afi? 2 2 



Muj 2 R 2 



E, 



where 



and 



V q ,x(R) = EP N -*>/*F q ,x{R) 



K ni = A + 



3N-6 



(25) 
(26) 



Noticing that the effective hyperradial Schrodinger equa- 
tion, Eq. (pM)) . is formally identical to the Schrodinger 
equation for the three-dimensional harmonic oscillator 
with angular momentum AT n i [20l [Hoj , the eigenenergies 
A ni , rc i, E^^N^NJ = £ ni _(AT t ,A^) - £ cm , and the 
corresponding eigenfunctions F q ,\(R) are readily written 
down, 



An 



with q = 0, 1, • • • and 



2? + AT n 



(27) 



N qKt .,R 



K ni + 1 



exp 



F 9 ,A(i?) 



(28) 



The normalization constant N q K ni is chosen such that 
|F 9 , A (A)| 2 dA=l, (29) 



TABLE I: Definition of the Jacobi vectors pi and the associ- 
ated reduced masses p.% used in our construction of the hyper- 
spherical harmonics $a,x f° r the two-component equal-mass 
atomic Fermi gas with (Nf, Aj,) = (2, 1) and (2, 2). 





Pl P2 f>3 


Mi 


/X2 




(2,1) 
(2,2) 


f - 

Tl - 


-f 2 £^-r 3 

- fi * 3 - ft ^ - 


m 
2 

m 

T 


2m 

3 
m 

2 


m 


leading to 










N qKni 










(30) 


= \ 


(2A- ni + l)!!0FAf" i+1/2) (O)« 


2Kni-t 
ho,M 


-3 ' 



where A 



(A' m + l/2) 



denotes the associated Laguerre poly- 




nomial. The harmonic oscillator length a^ ,Mi aho.M = 
\fh/ '(Mui), can be interpreted as being associated with 
an effective mass M particle that moves along the hy- 
perradial coordinate R. The quantity K n \ depends on A 
and can be thought of as an effective angular momen- 
tum quantum number; it should be noted, though, that 
the K nl are, in general, neither equal to the total an- 
gular momentum L tot nor equal to the relative angular 
momentum L re \ of the two-component Fermi gas. 

The explicit construction of the hyperspherical har- 
monics $a,x(^) requires the hyperangles Q to be speci- 
fied. The hyperangles f2 can be defined in many differ- 
ent ways and here we employ definitions that allow for a 
straightforward anti-symmetrization of the $a,x(^)- To 
this end, we introduce a set of mass scaled Jacobi vectors 
Hi, i = 1, • • • , N - 1, where 0, SF 



(31) 



Table Q] lists the Jacobi coordinates pi and the associ- 
ated reduced masses fa employed in this paper for the 
treatment of atomic two-component equal-mass three- 
and four-fermion systems. These Jacobi vectors have 
particularly convenient properties under the exchange of 
identical fermions. For the (A^, ATj.) = (2, 1) system, the 
Jacobi vector p\ changes sign under the exchange of the 
two spin-up fermions while pi remains unchanged. For 
the (Nf, Ni) — (2,2) system, the exchange of the two 
spin-up fermions leads to a sign change of p\ while p>2 
and ps remain unchanged, the exchange of the two spin- 
down fermions leads to a sign change of pi while p\ and 
P3 remain unchanged, and the simultaneous exchange of 
the two spin-up fermions and the two spin-down fermions 
leads to a sign change of p\ and p2 while p?, remains un- 
changed. These properties of the Jacobi vectors make 
the construction of properly anti-symmetrized hyper- 
spherical harmonics <3?a,x(^) f° r the three- and four- 
fermion systems comparatively simple. In terms of the 
mass-scaled Jacobi vectors Ui, the 3 A" — 4 hyperspher- 
ical angles SI are defined as Q = {ux/R, ■■■ ,un-x/R) 
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TABLE II: Characterization of the hyperspherical harmonics 
$a,x(^) for the ( N t, N i.) = (2,1) system with A < 5. In 
determining x> only hyperspherical harmonics that change 
sign under the exchange of the two spin-up atoms are counted. 



A 


X 


K ni 


L T c\ 


n rc i 


1 


3 


5/2 


1 


-l 


2 


5 


7/2 


2 


+i 


2 


3 


7/2 


1 


+i 


2 


1 


7/2 





+i 


3 


14 


9/2 


3 


-l 


3 


5 


9/2 


2 


-l 


3 


6 


9/2 


1 


-l 


4 


18 


11/2 


4 


+i 


4 


14 


11/2 


3 


+i 


4 


15 


11/2 


2 


+i 


4 


3 


11/2 


1 


+i 


4 


1 


11/2 





+i 


5 


33 


13/2 


5 


-l 


5 


18 


13/2 


4 


-l 


5 


28 


13/2 


3 


-l 


5 


10 


13/2 


2 


-l 


5 


9 


13/2 


1 


-l 



and the hyperradius R, Eq. (|18[) . can be rewritten as 

Following Avery [43 |. we construct a complete set of 
hyperspherical harmonics #a,x(^)j which are simulta- 
neous eigenfunctions of the operators A 2 , L 2 ol , L Te \ >z 
and LT rc i. Although the explicit functional forms of the 
<fr\ tX (£l) are needed in our perturbative treatment, we 
restrict ourselves here to summarizing the degeneracies 
of the non- interacting eigenfunctions for the (Nf, N±) — 
(2, 1) and (2, 2) systems (see Tables HJ and nj) . 

Tabic |TT] shows the degeneracies x an( i quantum num- 
bers of the hyperspherical harmonics for the (Nf, Nf) = 
(2, 1) system with A up to 5; in constructing Table ITT1 only 
hyperspherical harmonics that change sign under the ex- 
change of the two spin-up atoms have been counted. This 
symmetry constraint reduces the degeneracy of each A 
manifold tremendously. Equation (|2"3")) — applicable to a 
system without symmetry constraints — gives 1, 6, 20, 
and 50 for A = 0,1,2 and 3, while Table [TT1 shows that 
the degeneracies are reduced to 0, 3, 9 and 25. Table [TT1 
can be readily constructed by considering the angular 
momentum operators l\ and I2 associated with the Ja- 
cobi vectors pi and P2, and by taking into account that 
l\ and I2 couple to L rc \ [44|. Since pi is the Jacobi vec- 
tor that connects the two spin-up fermions, l\ can only 
take odd values; 1%, in contrast, is not restricted by sym- 
metry constraints, implying I2 — 0, 1, • • • . For a given 
A, the allowed (I1J2) combinations are determined by 
X = h + l 2 + 2p, where p = 0, 1, • • ■ @, [H|. Since the 
(^1,^2) = (0,0) combination is symmetry-forbidden, the 



TABLE III: Characterization of the hyperspherical harmonics 
$A, X (^) for the ( N t, N l) = ( 2 . 2 ) system with A < 5. In 



determining \, only hyperspherical harmonics that change 


sign 


under the exchange 


of the two 


spin-up atoms 


and the 


two 


spin-down atoms are 


counted. 






A 


X 


K ni 




Ilrel 


2 


5 


5 


2 


+1 


2 


3 


5 


1 


+1 


2 


1 


5 





+1 


3 


7 


6 


3 


-1 


3 


10 


6 


2 


-1 


3 


9 


6 


1 


-1 


3 


1 


6 





-1 


4 


27 


7 


4 


+1 


4 


28 


7 


3 


+1 


4 


35 


7 


2 


+1 


4 


12 


7 


1 


+1 


4 


3 


7 





+1 


5 


33 


8 


5 


-1 


5 


54 


8 


4 


-1 


5 


77 


8 


3 


-1 


5 


50 


8 


2 


-1 


5 


27 


8 


1 


-1 


5 


2 


8 





-1 



smallest allowed A value is 1. For A = 1, the only pos- 
sible (Zi,ia) combination is (1,0), resulting in L IC \ = 1, 
II rc i = (— l)'i+' 2 = —1 and a degeneracy of x — 3 (corre- 
sponding to three different projection quantum numbers 
Ml)- For A = 2, the only possible (I1J2) combination 
is (1,1), leading to L Te \ — 0,1,2 and n re i = +1. The 
degeneracy x is 9 (1, 3 and 5 states for L re i = 0, 1 and 
2, respectively). For A = 3, the allowed (Zi , Z2) combina- 
tions are (3,0), (1,2) and (1,0), leading to L Te \ = 3 (7 
states), L Ie \ — 3,2,1 (15 states) and L ie \ = 1 (3 states), 
respectively; thus, the degeneracy x is 25. Following this 
reasoning, the remaining entries in Table [TT] can be veri- 
fied. 

Table IIIII summarizes the degeneracies and quantum 
numbers for the (Nf,N±) = (2,2) system. Similarly to 
the three-fermion case, Table IIIII is constructed by re- 
alizing that the angular momentum quantum numbers 
h and I2 associated with the Jacobi vectors pi and f>2 
can only take odd values and that Z3, where Z3 denotes 
the angular momentum quantum number associated with 
the Jacobi vector p3, can take any value. For a given 
A, the allowed (h,l2,h) combinations are determine d by 
A = h + h + h + %P + 2q, where p, q = 0, 1, • • • 0, [H[. 
Since both l\ and 1% have to be odd, the smallest allowed 
A value is 2. In this case, l\ = I2 = 1 and Z3 = 0. The 
A = 2 manifold thus consists of nine states [Zi, I2 and 
l 3 can couple so that L re i = 2 (5 states), 1 (3 states) 
and (1 state)]. For A = 3, the only possibility is 
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TABLE IV: Coefficients c (1) for Fermi gas with 
(Nt,N±) = (2,1). The c (1) are defined through 
E (D =c (D(2 7 r)- 1 /2& tja M/ aho . 



-Eni.rol/ (fiw) 


Qni } xel 


L rc i 


n rc i 


c« 




4 


3 


1 


-l 


3 




5 


5 


2 


+i 


3/2 




5 


3 


1 


+i 







5 


1 





+i 


15/4 




6 


7 


3 


-l 


9/4 




6 


7 


3 


-l 







6 


5 


2 


-l 







6 


3 


1 


-l 


fg(13 + v 


/ 41) 


6 


3 


1 


-l 


£(13 -v 


/ 41) 


6 


3 


1 


-l 








{h,h,h) = (1, 1, 1), implying 27 states. The angular mo- 
menta corresponding to these 27 states can be obtained 
by first coupling l\ and I2 to an intermediate angular 
momentum vector with quantum number 2, 1 or 0, and 
then coupling the intermediate angular momentum vec- 
tor and I3 to obtain L IC \. The higher A manifolds are 
treated following the same scheme. 

Knowing the allowed A and \ values, the degeneracy 
5ni.roi of a given relative energy E n i iIC \(N-f, N^) of the non- 
interacting trapped Fermi gas can be easily determined 
using Eq. (|21)1) and Eq. (|2"T|) . These degeneracies are sum- 
marized in the second column of Tables IIVI and [V] for the 
three- and four-fermion systems, respectively. Alterna- 
tively [U the relative energy E nitIC \(N-f, ^Vj.) of the 
non-interacting three- and four-fermion systems can be 
written as E ni , lel {N t , JV4.) = E^VOj + l i + 3/2)fiw, 
where rij = 0, 1, • ■ ■ and where the allowed angular mo- 
mentum quantum numbers lj are determined by the sym- 
metry requirements (see above). Counting the possible 
combinations of lj and rij values and taking the (2lj + 1) 
degeneracy associated with each lj into account, gives the 
same results as those reported in the second column of 
Tables IIVI and [V] and also allows — using Eqs. (|21j|) and 
(|2T|) — for an independent determination of the A and x 
values given in the first two columns of Tables |TT] and Mil 

So far, we have discussed the hyperspherical frame- 
work for the non-interacting two-component Fermi gas. 
We now review the modifications needed when apply- 
ing this framework to the infinitely strongly-interacting 
unitary gas with zero-range two-body interactions. The 
zero-range two-body potential with infinite a^ aa ^ does not 
establish a meaningful length scale, leaving the oscillator 
length a no as the only length scale in the problem. Using 
scale invariance arguments, it has been shown 20] that a 
diverging s-wave scattering length a^ aa ^ implies that the 
wave function ijj nn (ri, • • • , fjv) at unitarity separates in 
the same way as that of the non-interacting system [see 
Eq. (|2"T])]. It follows that Eq. ([H]) applies not only to the 
non-interacting gas but also to the unitary gas if K n [ is 



TABLE V: Coefficients c (1) for Fermi gas with 
(iVfjiVj.) = (2,2). The c (1) are defined through 



E W =c (D( 27r ) 


- 1/2 huja (aa 


'/ a ho- 






R ■ „„i / (hin\ 


^?ni,rcl 




n,.„i 




13/2 


5 


2 


+ 1 


5 


13/2 


3 


1 


+1 


4 


13/2 


1 







13/2 


1 k lo 

10/ £ 


7 


Q 
O 


1 


7 /9 


15/2 


5 


2 


-i 


3 


15/2 


5 


2 


_i 


2 


15/2 


3 


1 


_i 


5 


15/2 


3 


1 


-i 


If 29 4- \/4T> 


15/2 


3 


1 


_i 




15/2 


1 





_i 





1 7/9 


q 


4 




g^' T V 'OJ 


1 7 /9 

1 l / £ 


■j 


4 

4; 




O j £ 


1 7 /9 
Li j £ 





A 


1 1 


1 f97 \ /7Q^ 


1 7/9 
1 1 1 £ 


7 





1 1 

1 1 


4 l y + V 1 'J 


1 7/9 

1 l / £ 


7 


Q 
O 


1 1 

1 1 


■3 



17/9 

1 l j £ 


7 


■) 

O 


1 1 


1 (Q *./Tf\ 


1 7/9 


7 


•i 

O 


4. 1 


n 


1 7/9 


r, 



2 






17/9 


r. 
o 


2 


11 




1 7/9 
if / £ 


r. 


2 


4. 1 


4.61321 


17/2 


5 


2 


T -1- 


3 217S3 


17/9 


r. 
o 


2 


4_1 


Q 91 ^40 


1 7/9 
x 1 1 £1 


r. 


2 


4. 1 


1 Q919Q 


1 7/9 
± 1 / £1 


r, 



2 


4. 1 


1 4^4^ 


1 7/9 


r. 



2 


4_1 


n 


1 7/9 


•■> 






1 1 

1 -L 




17/2 


3 


x 


4-1 


3 


17/2 


3 


1 


+1 


1.73167 


17/2 


3 


1 


+1 


0.512668 


17/2 


3 


1 


+1 





17/2 


1 





+1 


7.40848 


17/2 


1 





+1 


6.98138 


17/2 


1 





+1 


15/4 


17/2 


1 





+ 1 


2.89139 



replaced by -Kunit and if A is reinterpreted as the eigen- 
value of the hyperangular eigenequation that takes the 
two-body interactions into account. In the following, we 
use ifunit to denote the effective angular momentum of 
the unitary gas. Note that K unit depends on the eigen- 
value of the hyperangular eigenequation, i.e., there exists 
a -Kunit for each channel function $a,x(^)' f° r notational 
convenience, we do not explicitly indicate the dependence 
of -Kunit on the hyperangular quantum numbers. 

The coefficients -Kunit have beeen obtained for all states 
of the three-fermion system [32j (see also Refs. (62T - 
Hi| for earlier work) and for the lowest 20 states with 
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(L re i,n re i) = (0, +1) of the four-fermion system 65[ by 



solving the hyperangular Schrodinger equation that in- 
cludes the two-body interactions (see also Ref. [42j ) . The 
relative eigenenergies -E U nit,rci of the unitary gas are, sim- 
ilarly to the non-interacting case, given by [20| 



E 



unit.rcl 



(2q + K unit + 3/2)fiw, 



(32) 



and Eqs. ([28|) - ([30|) remain valid if K n \ is replaced by -Kunit 
(and if A is reinterpreted as discussed above). In Sec. lIII| 
we determine a number of K unlt coefficients by solving 
the full relative Schrodinger equation of the four-fermion 
system for various L ro i > and by then comparing the 
resulting energy with the right hand side of Eq. (|3"2"1) . 
Lastly, we note that Eq. (|32j) implies that the excitation 
spectrum of the trapped unitary two-component Fermi 
gas contains ladders of excitation frequencies that are 
integer multiples of 2fiw, independent of the actual values 

Of if unit 0, El, Hil. 



C. Numerical treatment of the three-fermion 
system: Lippmann-Schwinger equation 

The trapped three-fermion problem with zero-range in- 
teractions and arbitrary s-wave scattering length a( aa ' 
has been solved using a number of different semi- 
analytical and numerical approaches (33l - [35| . Here, 
we replace the regularized zero-range pseudo-potential 
V^"gg(r*), which describes the interactions between atoms 
with opposite spins, by the corresponding Bethe-Peierls 
boundary condition and employ an approach developed 
by Kestner and Duan [33[ that is based on the Lippmann- 
Schwinger equation. This approach reduces the three- 
body problem to solving a set of coupled equations [33| , 



2r(-^-) 
1 (-»*-% 



B 

£< 



(act) » 
-,(aa) 



(33) 



value (a (Qa) /a„o Q) 



c = (ci 
In Eq. 



,c.b), and the eigen- 
, the Vj denote non- 



for the eigenvector c, 

i /t , - 1 . In Eq. (J33J), 
integer quantum numbers that depend on E ie \(2, 1) and 
the djk dimensionless matrix elements. Their definitions 
are given in Refs. [33l. lo6j. 

In the B — > oo limit, Eq. (f3"3"f gives the exact three- 
fermion energy spectrum. For each i? re i(2,l), there ex- 
ist multiple c and a (aa) that solve Eq. ([55]). Thus, 



Eq. (I33p can be interpreted as a matrix equation with 
eigenvector matrix c = (ci, ■ • • ,cb) and eigenvalue vector 
((a (aa) )i, • • • , (a (aa) ) B ). The solutions obtained by solv- 
ing Eq. (|33l) belong to three-fermion states with natural 
parity. For the three-fermion system, unnatural parity 
states are not affected by the s-wave zero-range inter- 
actions and coincide with those of the non-interacting 
system. For each L re i, we solve the matrix problem for 
different B, B < 50. For positive E Te i(2, 1), our results 
presented in Sec. IIII Al are obtained using B = 50. For 
negative -E rc i(2, 1), we use somewhat smaller B values; we 



have checked through extrapolation to the B — > oo limit 
that the three-fermion energies obtained in this manner 
are highly accurate. We find, e.g., that the eigenenergies 
at unitarity obtained by the numerical approach based 
on the Lippmann-Schwinger equation f33J agree to better 
than 0.01% with those obtained by solving the transcen- 
dental equation derived by Werner and Castin [32|. 



D. Numerical treatment of the four-fermion 
system: Stochastic variational approach 

To determine the energy spectrum of two-component 
Fermi gases with (Nf, N±) — (2, 2) under spherically sym- 
metric harmonic confinement, we employ the stochastic 
variational approach [35T - l4lj ]. Our implementation sep- 
arates off the center-of-mass degrees of freedom R cm , 
defines a set of N — 1 Jacobi coordinates x, x = 
(pi, ■ ■ ■ , /Ojv— i), and expands the relative wave function 
ip le \(x) in terms of the basis functions <Pk{i*), 



B 

E 

fe=i 



A[c k ifik (x)} 



(34) 



where the anti-symmetrization operator A can be written 
as A = 1 - P12 - P 34 + P12P34 for the (JV t , N±) = (2, 2) 
system. In Eq. (|34p . the c k denote expansion coeffi- 
cients. We parametrize the two-body potential Vtb by 
a spherically symmetric attractive Gaussian with depth 
Vq (Vq > 0) and range r , 



V G (r) 



-Vbexp 



\/2> 



(35) 



this interaction potential is convenient since the matrix 
elements (ifj\VQ\(pk) are — for the <fk employed in this 
work — known analytically [38[ (see below for the defi- 
nition of the ipk)- For a given range r , we adjust the 
depth Vo such that Vq reproduces the desired s-wave 
scattering length cS aa \ For negative (positive) a^ aa \ we 
restrict ourselves to parameter combinations for which 
Voir) supports no (one) s-wave free-space bound state. 
For each a^ aa ^ , we consider a number of different ranges 
ro, ro <C a n o, and extrapolate to the tq — > limit (see 
Sec.|InB]for details). 

We employ two different classes of non-orthogonal ba- 
sis functions tpk- For both classes, the Hamiltonian ma- 
trix elements and overlap matrix elements are known 
analytically (38|, thus reducing the problem of finding 
the eigenenergies and eigenfunctions to diagonalizing a 
generalized eigenvalue problem. The first class of basis 
functions Lp k (x) has well defined angular momentum L [e \ 
and natural parity, while the latter class has neither well 
defined angular momentum L rc \ nor well defined parity 

n re i. 

To describe natural parity states with well defined 
angular momentum L rc \ and corresponding projection 
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quantum number Ml, we employ the following basis 
functions [Hj], 

<p k (g) = \^\ L ^Y LiclML (v^)exp (^- l -x T A^x^ ,(36) 



where 



N-l 



,(k) 



(37) 



Here, the Uj, j 



l,--- ,JV - 1, define a (TV - 1)- 
dimensional parameter vector that determines how the 
angular momentum L IC \ is distributed among the (N — 1) 
Jacobi vectors pj. In Eq. (1551) . -A^ denotes a (TV — 
1) x (AT — l)-dimensional symmetric matrix, which is de- 
scribed by N(N — l)/2 independent parameters. To get 
a physical interpretation of these parameters, we rewrite 
the exponent on the right hand side of Eq. (|3f)|) in terms 
of a sum over the square of interparticle distances ry and 

N(N- l)/2 widths d['' 



O) 



AT-l JV-1 

•) ^ "ft 

i=l j=l 




(38) 



The explicit relationship between the parameter matrix 
A {k) and the widths d\f (i < j) can be determined 
by expressing the interparticle distance vectors ?ij in 
terms of the Jacobi vectors x [38}. Equation ([35]) il- 

(k) 

lustrates that the d)-' determine the widths of Gaus- 
sian functions in the interparticle distance coordinates. 

(k) 

In our calculations, we choose a set of widths d\j for 
each basis function and construct the matrix A^ k ' from 



(k) 

these. The widths d\j themselves are — guided by physi- 
cal arguments — determined semi-stochastically following 
the schemes discussed in Refs. [H, [H, HjJ. For the 
(N^,Ni) = (2,2) system with small ro and small posi- 
tive a,( aa \ e .g.. three-body and four-body bound states 
are absent [U [H[, implying that at most two of the 
widths d[ k ^ , d[ k ^ , d^ and d^ (but not d[ k ^ and d^ si- 
multaneously or d^ and simultaneously) should be 
of the order of the two-body range ro for a given k. We 
use the basis functions given in Eq. (J36J) to determine the 
eigenenergies of states with vanishing and finite L re \ and 
natural parity, i.e., n rc i = (— l) Lrcl . 

To describe states with unnatural parity, we employ 
basis functions ifk that are neither cigenfunctions of the 
angular momentum operator L re \ nor the parity operator 

n re i [H, 



ip k (x) = exp 



~x T A^x+(^) T x 



(39) 



Here, the quantity s^ fc ' consists of N — 1 three- 
dimensional parameter vectors and (s^ kS> ) T x is just the 
dot product between two 3(N — 1) dimensional vectors. 
The 3(7V — 1) parameters of are, together with the 
N(N — l)/2 parameters of the matrix A^ k \ optimized 
semi-stochastically [38]. Since the basis functions de- 
fined in Eq. (|3"9"|) are neither eigenfunctions of L re i, Ml 
nor n ro i , their use allows for the determination of the en- 
tire energy spectrum at once. In Sec. IIII Bl we employ 
the basis functions given in Eq. (|39[) to determine the 
energetically lowest lying unnatural parity state of the 
four-fermion system with negative a^ aa \ 

Following the schemes outlined, the determination of 
the four-fermion energies corresponding to unnatural 
parity states is significantly less numerically efficient than 
that of natural parity states. This is, of course, not 
sursprising since only a "fraction" of the basis functions 
given in Eq. (|3"9")l contributes to describing states with the 
desired angular momentum, projection quantum number 
and parity. 



III. RESULTS 

This section summarizes the energetics of the three- 
and four-particle equal-mass Fermi gas. 



A. Three-fermion system 

Figures [^a) and (b) show the eigenenergies E re \(2, 1) 
for states with natural parity and L rc i = and 1 as a 
function of the inverse scattering length l/a^ aa \ The 
symbols show the solutions to the coupled equations, 
Eq. (|33|) . while the solid and dashed lines are obtained 
from our perturbative treatments of the atomic Fermi 
gas and the effective atom plus dimcr model, respectively. 
Note that the perturbative treatment of the atomic Fermi 
gas (solid lines) describes the energy levels corresponding 
to gas-like states for negative as well as for positive a^ aa ' 
(|a( aQ )| small). In the following, we highlight selected 
characteristics of the three-fermion energy spectrum. 

We first consider the weakly-interacting attractive 
Fermi gas. In the non-interacting limit, a^ aa ^ — > CP, 
the ground state has an energy of -E n i,rei = 4fiw and is 
characterized by L rc i = 1 and n rc i = —1 (see Fig. [2] 
and Table IIV|) . For the next family of energies with 
£"ni.rci = 5hu>, we have two natural parity states with 
L le i = and L rc \ = 2, respectively, and one unnatural 
parity state with L rc \ — 1. The fact that the lowest non- 
interacting L IC \ = state has a higher energy than the 
lowest non-interacting L rc i = 1 state can be understood 
intuitively by realizing that the two like atoms cannot 
both occupy the lowest single particle state. Within the 
hyperspherical description, this implies that the L rc \ = 
state with A = and q — is symmetry-forbidden (see 
Sec. Ill Bp and that the first symmetry- allowed L re i = 



10 




FIG. 2: (Color online) Three-fermion energies for (a) 
(Z/rel,n re i) = (0, +1) and (b) (L Te i, U le \) = (1,-1) as a func- 
tion of the inverse s-wave scattering length a\ lo /a^ aa \ Sym- 
bols show the essentially exact zero-range energies obtained 
by solving Eq. (|33[1 . Solid lines show the energies obtained by 
treating the non-interacting atomic Fermi gas perturbatively 
for negative and positive a' aa ' while dashed lines show the 
energies obtained by treating the effective atom plus dimer 
system perturbatively. 



state, which has A = 2 and q — 0, lies 2huj higher in 
energy than the symmetry-forbidden L rc \ = state. 

The coefficients c^ 1 ' that determine the perturbative 
energy shifts E^ of the atomic (iV-f, N±) = (2,1) sys- 
tem are calculated semi- analytically following the scheme 
outlined in Sees. Ill Al and III Bt and reported in the last 
column of Table IIVI for the first three energy families. 
As already mentioned, the unnatural parity states of 
the three-fermion system are uneffected by the interac- 
tions [2(J HH HH . This behavior is specific to zero-range 
s-wave interactions since a finite-range potential allows, 
in general, for an energy shift due to p-wave or other 
higher partial wave interactions. To illustrate the valid- 
ity regime of the perturbative expressions, Figs. EJa)-(d) 
show the small |a^ aa ^| region, a' aa ) < 0, of the three- 



region, 

body energy spectrum as a function of 



(aa)| 



aho 



for 



the energies around the first four energy families with 
E Tcl (2, 1) « 4fuv to E Tel (2, 1) « Ihjj. As in Fig.d the ex- 
act energies are shown by symbols while the perturbative 
energies corresponding to natural parity states are shown 
by solid lines. As expected, the perturbative treatment 
reproduces the exact energies extremely well for small 




0.05 0.1 0.15 0.2 

, (aa). 



/a 



ho 



FIG. 3: (Color online) Three-fermion energies as a function 
of the absolute value of the s-wave scattering length |a^ aa ^| 
for small |a (aa) |, a (aa) < 0, and (a) £ re i(2, 1) « 4Hlj, (b) 
£ r el(2,l) « bhio, (c) £W(2,1) w 6fiw, and (d) E re i(2, 1) w 
7fiw. Squares, pluses, diamonds, crosses and circles show 
the essentially exact zero-range energies obtained by solving 
Eq. (|33[) for L rc i = — 4. For comparison, solid lines show the 
perturbative results for natural parity states. The solid lines 
correspond to (a) (Z/ re i,n rc i) = (1,-1); and from bottom to 
top to (b) (L rel ,II rcl ) = (0,+l) and (2,+l); (c) (L rcl ,II rcl ) = 
(1,-1), (3,-1), and (1,-1); and (d) (L re i,n re i) = (0,+l), 
(2,+l), (0,+l), (4,+l), and (2,+l). 
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FIG. 4: Frequency of the coefficients c- , which character- 
ize the energy shifts E^ 1 ' of the weakly-interacting atomic 
Fermi gas, for the three- and four-fermion systems. Pan- 
els (a) and (b) show the distribution of the c' 1 ' coefficients 
for E le i(2, 1) « 7huj (fourth energy manifold of the three- 
fermion system) and for E rc \(2, 1) w 8huu (fifth energy man- 
ifold of the three- fermion system). Panels (c) and (d) show 
the distribution of the coefficients for E xe i(2,2) « 19ftu/2 
(fourth energy manifold of the four-fermion system) and for 
E Te \(2, 2) ss 21fiw/2 (fifth energy manifold of the four-fermion 
system). Note the log scale of the vertical axis. 

|a( aa ) | / aho and provides a semi-quantitatively correct de- 
scription up to |et( aa '| 0.5aho (see also Fig. [21 which 
shows the perturbative energies up to aho/|c^ aa '| = 2 or 
\a<- a ^\ =0.5a ho ). 

Tabic ITVl shows that the coefficients cover a wide 
range of values. Within each energy family, the state 
shifted most strongly by the interactions is a natural 
parity state with the smallest allowed angular momen- 
tum L Te \. To illustrate the increasing density of states 
and the spread of the energy levels around the non- 
interacting degenerate energy manifold, Figs. Ufa) and 
(b) show the frequency with which the coefficients c' 1 ) 
occur for E re \(2, 1) w 7huj (fourth energy manifold) and 
E xe \(2, 1) ~ 8hus (fifth energy manifold), respectively. In 
making this plot, the 2L rc \ + 1 degeneracy of the energy 
levels has been taken into account. Since the unnatural 
parity states are not affected by the zero-range interac- 
tions, the distribution of the coefficients shows a large 
amplitude for the = bin. Figures@|a) and (b) show 
that the spread of the coefficients increases slightly as 
the energy manifold increases. The primary characteris- 
tic of the distributions of the cW coefficients is, however, 
that the amplitude increases with increasing energy. 

We now consider the weakly- repulsive regime, i.e., the 
regime where a^ aa ^/aho *C 1 and a( aa ) > 0. In this 
regime, Fig. [2] shows two families of energy levels: (i) 
energy levels with positive energy and (ii) energy levels 
with negative energy. The positive energy branches cor- 
respond to states that describe a gas of atoms; we refer 
to these states as the "gas-like state" family. The ener- 
gies of this family are, in the a^ aa ^ —> + limit, well de- 



scribed by treating the atomic Fermi gas perturbatively 
(see solid lines in Fig. [5J . The negative energy branches 
correspond to states that can be thought of as consisting 
of a bound diatomic molecule and a spare atom; we refer 
to these states as "dimer plus atom" family. In agree- 
ment with the literature (see, e.g., Ref. [561]). Fig.[5]shows 
that the formation of bound triatomic molecules is pro- 
hibited by the Pauli exclusion principle or the so-called 
Pauli pressure. For small and positive a' oa ', the pertur- 
bative energy shifts for the energy levels with L xc \ = 
and II ro i = +1 [dashed lines in Fig. Ufa)] are calculated 
using Eq. (|17[) with k = ad. The perturbative approach, 
applied to the effective model Hamiltonian H , predicts 
no energy shift for states with L Te \ > 0. This follows 
directly from the fact that we parametrized the effec- 
tive atom-dimer interaction through a zero-range s-wave 
potential. Thus, the "bending" of the dashed lines in 
Fig. [2fb) for L IC \ = 1 (and in general, L rc i > 0) and 
positive cS aa } is solely due to the internal energy of the 
dimer and not due to the effective atom-dimer interac- 
tion. We find that the perturbative treatment provides 
a qualitatively correct description up to « 0.5a ho . 

The effective model Hamiltonian H eS also provides an 
intuitive picture for why the lowest L rc \ — state has a 
lower energy than the lowest L re i = 1 state as a^ aa ^ — » + . 
In the a( aaS> — > + limit, the diatomic molecule has van- 
ishing angular momentum and the angular momentum 
L re ] must be carried by the atom-dimer distance vector. 
Thus, the energy of the lowest state with L Ie ] = 1 lies 
approximately Hui above the lowest state with L Te \ = 0, 
the energy of the lowest state with L rc \ = 2 lies approx- 
imately hx> above the lowest state with L Te \ — 1, and so 
on. The parity inversion of the energetically lowest lying 
state (L Te i — 1 and II rc i = —1 in the a^ aa " > — > CP limit, 
and L Ie i = and n re i = +1 in the a (aa) -> 0+ limit) oc- 
curs at aho/a^ aa ' ~ 1 and has already been pointed out 
by a number of works (33rl35j . 

Motivated by the fact that the energy spectrum of the 
three-fermion system can be described by an effective 
atom plus dimer model in the a^ aa ' — > + limit, symbols 
in Figs. Efa)-(d) show the quantity E le \(2, 1) — Bdimer 
as a function of the inverse atom-atom scattering length 
l/ a O a ) f or £ rel = — 3. Here, E'dimcr denotes the low- 
est eigenenergy of the trapped s-wave interacting atom- 
atom system, i.e., the lowest eigenenergy of Eq. (| 15[) with 
k = aa. The quantity E ve \(2, 1) — Sdimor has been investi- 
gated previously [H, HH and has been termed the univer- 
sal energy crossover curve in Ref. (35|. Figures Efb)-(d) 
show the existence of a family of states for L IC \ > whose 
scaled energies are nearly independent of the s-wave scat- 
tering length a,( aa \ These scaled energies are approxi- 
mately given by (2n c g + L re \ + 3/2)hu [see dashed lines 
in Figs. Efb)-(d)]. The fairly good agreement between 
the symbols and the dashed lines reflects the fact that a 
subset of the three-fermion energies can be described to a 
fairly good approximation by treating the three-fermion 
system as consisting of a bound trapped s-wave dimer 
plus a non-interacting spare atom. Figures[5Jb)-(d) show 
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FIG. 5: (Color online) Symbols show the scaled energies 
E Te \(2, 1) — -Edimer as a function of the inverse scattering length 
a ho /a^ for (a) (L rel ,Il rcl ) = (0, +1), (b) (L rcl ,Il rcl ) = 
(1,-1), (c) (L re i,n re i) = (2,+l), and (d) (L rc i,n rc i) = 
(3, —1). For comparison, dashed lines show the energies pre- 
dicted by the effective atom-dimer model, which provides a 
semi-quantitative description of the states belonging to the 
"atom plus dimer" family (see text for details). 



that this effective dimer-atom description improves with 
increasing L le i- 

The fact that a subset of states of the three-particle 
spectrum is reasonably well described by the effective 
dimer plus atom model suggests that the three-particle 
energy spectrum can be described in terms of avoided 



crossings between "atom plus dimer" states and "gas- 
like" states. An interpretation along this line has been 
quantified by von Stecher [67} who applied a diabatiza- 
tion scheme. Here, we do not follow the diabatization 
scheme but instead offer a qualitative discussion of the 
natural parity three-fermion spectrum with L rc \ > 0. We 
make four observations: (i) A sequence of states has an 
energy E Ie i(2, 1) of approximately (2 + L re i + 2n)fkv at 
unitarity, an energy £ , re i(2, 1) of approximately (3+£ re i+ 
2n)huj in the a^ aa ^ — > 0~ limit, and an energy £" re i(2, 1) 
of approximately -Edimcr + (3/2 + L re i + 2n)hw in the 
a (aa) _^ q+ ij m }t j see pjg fjjb); the states discussed here 
are those with approximately constant E re \(2, 1) — -Edimer, 
see Figs. |S[b)-(d)]. (ii) For each L ro i, the degeneracy 
of the energy families with -E n i,rei = (3 + L rc \ + 2n)huj, 
n = 0, 1, • • • , increases by one (or 2L le \ + 1 if the de- 
generacy of the different Ml values is accounted for ex- 
plicitly) as n increases by one (see Fig. [2J Fig. [5] and Ta- 
blelIV[). (iii) The energy spectrum corresponding to "gas- 
like" states has to be identical in the limits cS aa ^ — > 0~ 
ando'""' 0+. (iv) It can be easily checked that (i)-(iii) 
are consistent with the fact that the energy K re i(2, 1) of 
all but one level of each non-interacting manifold with a 
given L rc \ decreases by hw when going from a^ aa ^ — > 0~ 
to a( QQ ) -> oo and by another hjj when going from 
a {aa) _^ oo to a^ aa ) — » + . The dropping of the ener- 
gies by 2huj as a^ aa ^ changes from 0~ through ±oo to + 
is similar to the 2hw dropping of the excited state s-wave 
energies of the two-particle system (see dashed lines in 
Fig. [I]) and can be interpreted as one pair (consisting 
of a spin-up atom and a spin-down atom) feeling the s- 
wave interaction while the other spin-up atom carrying 
the angular momentum. 

Interestingly, the observations described in the previ- 
ous paragraph allow for an approximate determination 
of the .Kunit coefficients [see Eq. ([32)) ]. Observation (i) 
implies that the lowest -K un it coefficient for a given L Te \ 
is approximately given by -K umt; modci = L xe \ + 1/2 in the 
large L rc i limit. As discussed in Sec. Ill Bl each -K un i t coef- 
ficient determines the starting point of a ladder of energy 
levels, which are spaced by 2hjj and which are associated 
with an increasing number of nodes along the hyperradial 
coordinate R. It is evident from Figs.[HJb)-(d) that these 
states, which are characterized by the same hyperangu- 
lar quantum number but different hyperradial quantum 
numbers q, transform into atom plus dimer states in the 
a {aa) Q+ limit, which are characterized by the effec- 
tive orbital angular momentum quantum number l e g = 
and different radial quantum numbers n c g. Interpreting 
the atom-dimer distance coordinate as the hyperradial 
coordinate in the a t - aa ^ — > + limit, the identification 
q O n e ff suggests itself. Using observation (iv), the re- 
maining -K^it coefficients at unitarity are approximately 
given by K" U mt,model = L xe \+l/2+2n for L IC \ > and nat- 
ural parity states [for each n = 1, 2, • • • , the q quantum 
number in Eq. (l32|) takes the values q — 0, 1, • • •]. Fig- 
ure [6ja) shows that the difference between K un a (sym- 
bols) and -Kunit, model (dotted lines) decreases as L ro \ in- 
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FIG. 6: (Color online) Symbols show the coefficients K un it 
that correspond to natural parity states for (a) the three- 
fermion system and (b) the four-fermion system as a function 
of the angular momentum quantum number L te \. Dotted lines 
show the coefficients .Kunit, model (see text for details). 



creases. The difference between -K un j t and if unit, model has 
previously been quantified by Werner and Castin within 
a semi-classical WKB framework [32] . 

The L rc \ = spectrum is different from the L re \ > 
spectra for two reasons. First, the effective atom-dimer 
system is described by an effective atom-dimer s-wave 
scattering length a^ ad \ which leads to a decrease of ap- 
proximately fko of the energy levels belonging to the atom 
plus dimer family as a( aa ' changes from oo to + . Sec- 
ond, the symmetry constraint in the a^ aa ^ — > 0~ limit 
pushes the energy of the lowest L IC \ — state up by huj 
compared to that with L rc i = 1. Dashed lines in Fig.[5ja) 
show the eigenenergies of the effective atom plus dimer 
model, i.e., the eigenenergies of Eq. (TTSl) for k = ad with 
a ( ad ) given by Eq. (fT3|) . It can be seen that the agreement 
between the energies of this effective model and a subset 
of the full three-fermion energies is good for small posi- 
tive a^ 00 ) and qualitatively correct throughout the entire 
crossover regime. Using the effective atom plus dimer 
model, the energy at unitarity of a subset of states is ap- 
proximately given by -E re i(2, 1) « -Edimer + (2q + 5/2)fiw, 
implying -Kunit. model — 3/2. For comparison, the exact 
value is -Kunit = 1.666 [32J. The other -K un i t coefficients 
for L rc \ = can be estimated by using that the energy of 



to a (aa ) -> 0+, implying £ U nit,rci(2, 1) ~ {5/2+2q+2n)huj 
with n = 1, 2, • • • and q = 0, 1, • • • . This estimate yields 
-Ku„it,modei = 5/2 + 2n for n = 1, 2, ■ ■ • . Figure EJa) 
shows that the -K un it, model (dotted lines) reproduce the 
exact -K un it coefficients (symbols) fairly well. 



B. Four-fermion system 

This section discusses the energy spectrum of the four- 
fermion system throughout the BCS-BEC crossover. We 
primarily focus on the energies corresponding to nat- 
ural parity states but also consider those correspond- 
ing to unnatural states. To determine the energy spec- 
trum corresponding to states with natural parity, we use 
the stochastic variational approach with the basis func- 
tions given in Eq. (1361) ; our basis set optimization ei- 
ther treats one state at a time or a subset of states si- 
multaneously. For a given atom-atom scattering length 



, we determine the energies for various ranges ro, 



ro = 0.01a no — 0.09aho, of the two-body Gaussian inter- 
action potential and extrapolate the finite-range energies 
to ro — > 0. For negative scattering lengths a i ~ aa \ we find 
that the four-fermion energies depend linearly on ro for 
all I/ re i considered. In this regime, we typically calculate 
the energies for three different r and then determine the 
ro — > energies by performing a linear fit. For positive 
scattering lengths cS aa \ the energetically low-lying part 
of the spectrum is dominated by the "internal" energy of 
the dimer(s) formed. As discussed in more detail below, 
the lowest energy family for even L IC \ can be described 
by an effective two-boson model while the lowest energy 
family for odd L IC \ consists of states that can be thought 
of as consisting of a dimer and two atoms as cS aa ^ — > + 
(see Sec. |TT] and below). Correspondingly, for positive 
a (aa) anc j even £ rclj we subtract twice the dimer binding 
energy from the four-fermion energies for each ro and ex- 
trapolate the scaled four-fermion energies to the ro ^ 
limit. We typically consider five different ro and extract 
the scaled zero-range energies by performing a quadratic 
fit to the scaled finite-range four-fermion energies. The 
zero-range four-fermion energies themselves are then ob- 
tained by adding twice the zero-range dimer energy, i.e., 
2-Kdi mcr . For odd L re i, we subtract (and later add) the 
dimer energy as opposed to twice the dimer energy but 
proceed analogously otherwise. 

The energies of a large number of energy levels cor- 
responding to natural parity states are reported in the 
auxiliary materials (68| for atom-atom scattering lengths 



ranging from over oo to 0.2a no for L IC \ < 4. To 



one subset of states drops by fku in going from a 



(aa) 



OO 



the best of our knowledge, these are the first comprehen- 
sive benchmark results for the four-fermion system with 
finite angular momentum throughout the crossover re- 
gion. Tables IVII and I VIII summarize the energies of the 
four-fermion system at unitarity for natural parity states 
with L rc \ = to L rc \ = 8 as well as for one unnatural 
parity state. For L YC \ > 0, these are the first results at 
unitarity. We estimate that our extrapolated zero-range 
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TABLE VI: Extrapolated zero-range energies _E un it,rei(2, 2) for 
natural parity states with L re j < 4 [£7 un it,rci(2, 2) < 8.57k*;]. 
The uncertainty of the energies is estimated to be in the 
last digit reported. For comparison, the lowest energy of the 
(L reI ,n re i) = (1,+1) state is £ U nit,rei = 5.088(20)&j. 
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TABLE VII: Lowest two extrapolated zero-range energies 
-E U nit,rci(2, 2) for the four-fermion system with natural parity 
and L rc i = 5 — 8 at unitarity. The uncertainty of the extrap- 
olated zero-range energies is estimated to be in the last or 
second last digit reported. 
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energies for the energetically lowest-lying L rc \ = state is 
accurate to better than 0.1% for most scattering lengths 
a (aa) ^ mcm ding infinitely large a^ aa \ Near the avoided 
crossings around a^ aa ^ ~ aho (see, e.g., Fig.Qjll), however, 
the accuracy decreases by up to an order of magnitude. 
Generally speaking, we find that the accuracy of the ex- 
trapolated zero-range energies also decreases for energet- 
ically higher lying states and for states with larger L Te \. 
The eigenenergies reported in Tables I VII and IVIII are la- 
beled by the hyperradial quantum number q. Following 
Eq. p2[) . we identify this quantum number by looking 
for 2fvjj spacings between energy pairs. Inspection of Ta- 
ble |VI] shows that the energies with q > lie, within our 
numerical accuracy, 2hiu above an energy with q — 1 (see 
the nearly identical K un a coefficients in the fourth col- 
umn of Tables I VII and fVIIf) . Figure [(Jb) shows the K unit 
coefficients corresponding to natural parity states of the 
four-fermion system as a function of L IC \. Compared to 
the three-fermion system, the four-fermion system ex- 
hibits a notably denser energy spectrum at unitarity [see 
Fig. Efa) and Mb)]. 

Lines in Fig. [7] show the extrapolated zero-range ener- 
gies for the four-fermion system corresponding to natu- 
ral parity states with L rc \ = to L rc \ = 4 as a function 
of the inverse s-wave scattering length l/a^ aa ' for nega- 
tive a^ aa \ In the a^ aa ' — > CP limit, the three energeti- 
cally lowest-lying four-fermion energy manifolds around 
£ni,rci(2,2) = 13ftw/2, and 177kj/2 consist of two, 

four and 15 states, respectively (here, the 2L IC \ + 1 de- 
generacy due to the Ml quantum number is not included 
in counting the states; see also Table N} . For compar- 
ison, pluses in Fig. [7] show the energetically lowest ly- 
ing unnatural parity state with (L rc i,II rc i) = (1,+1); 
in this case, the energies are calculated for a Gaussian 
two-body potential with small but finite r and have not 
been extrapolated to the ro — > limit. Figure [JJ shows 
that the three energy manifolds remain distinguishable 
up to a\ lo /cS aa } w —2 but start to overlap notably in the 
strongly-interacting regime. 

We now discuss the weakly-attractive regime of the 
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FIG. 7: (Color online) Four-fermion energies -E re i(2,2) as a 
function of the inverse s-wave scattering length aho/a' aa ' with 
l/a (aa) < 0. For a (aa) -» CT, all energy levels correspond- 
ing to natural parity states around -E n j ir ci = 13fiu/2, 15ftu/2 
and YThw/2 are shown. Solid, dotted, dashed, dash-dotted 
and dash-dash-dotted lines show the zero-range energies cor- 
responding to natural parity states with L TC \ = to L r ci = 4. 
Pluses show the energy of the energetically lowest-lying un- 
natural parity state with L rc \ = 1. 




FIG. 8: (Color online) Four-fermion energies _£/ rc i(2,2) as 
a function of the absolute value of the s-wave scattering 
length \a {aa) \ for small |a (aa) |, a (aa) < 0, and £ rcl (2,2) w 
13huj/2. Squares, diamond and pluses show the numerically 
determined four-fermion energies for (L rc i,II ro i) = (0,-1-1), 
(L rc i,n rc i) = (2,+l) and (L Ieh n re i) = (1,+1), respectively. 
Solid and dashed lines show the perturbative results for nat- 
ural and unnatural parity states, respectively. 



four-fermion energy spectrum in more detail. The coeffi- 
cients c^ 1 - 1 that determine the perturbative energy shifts 
E^ are summarized in the last column of Table [V] for 
the first three energy manifolds with £ni. r oi = l3frio/2, 
Xhhbjjl and 17haj/2. Figure [S] compares the four-fermion 
spectrum near E rc \(2,2) w 13/2huj calculated by the 
stochastic variational approach [squares, diamonds and 




FIG. 9: (Color online) Four-fermion energies -E/ re i(2, 2) as a 
function of the absolute value of the s-wave scattering length 
|a (aa) | for small \a (aa) \ (a (aa) < 0), £ rc i(2,2) « 17fia;/2 and 
(L rc i,II ro i) = (2, +1). Diamonds show the numerically deter- 
mined four-fermion energies while solid lines show the pertur- 
bative results. Note that the fourth and fifth state (counted 
from the bottom) are nearly degenerate for the scattering 
length range depicted. 

pluses show the energy levels corresponding to states with 
(L ro i,n ro i) = (0,+l), (2,-1-1) and (1,+1), respectively] 
with that calculated perturbatively (solid and dashed 
lines show the energies corresponding to states with natu- 
ral and unnatural parity, respectively). As expected, the 
agreement is excellent for small |a^ aa - ) | and worsens with 
increasing |o^°°^|. For small |a^°°^|, the energy level with 
unnatural parity is affected less strongly by the two-body 
interactions than the energy levels with natural parity. 
Inspection of Table [Vl shows that this is a general trend, 
i.e., within a given manifold the energy level shifted most 
strongly is that corresponding to the natural parity state 
with the smallest allowed angular momentum. 

To illustrate the behavior of the energy spectrum for 
a higher energy manifold, Fig. [9] shows the energies cor- 
responding to the eight states with (i re i,n re i) = (2, +1) 
around E Te \ rs Y!hw/2. Again, the agreement between 
the numerically determined energies (diamonds) and the 
perturbatively determined energies (solid lines) is ex- 
cellent for small |a ( - aa -'|. Interestingly, within the per- 
turbative treatment, one of the (L rc i,n rc i) = (2, +1) 
states is not affected by the s-wave interactions, imply- 
ing that the wave function vanishes whenever two un- 
like fermions approach each other closely. It turns out 
that the perturbative result in this case is exact, i.e., 
there exists a (L re i, n rc i) = (2, +1) state with energy 
£ rc i(2,2) = nhw/2 for all scattering lengths a {aa) . 

FiguresHfc) and (d) show the distributions of the c' 1 - 1 
coefficients for the fourth and fifth energy manifolds with 
# re l(2,2) rj 19ftw/2 and E re i(2, 2) « 21ftw/2, respec- 
tively. Compared to the three-particle case [Figs. Ufa) 
and (b)], the degeneracies increase more rapidly as can 
be seen by the higher frequency with which the c*- 1 ' co- 
efficients occur. Furthermore, the c^ 1 ' = bin no longer 
dominates the distribution since both natural and unnat- 
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FIG. 10: (Color online) Scaled four-fermion energies 
E Te i(2, 2) — 2E'dimcr as a function of ai 10 /a' aa ' for (a) 
(ird,n rc i) = (0,+l), (b) (L re i,n re i) = (2,+l) and (c) 
(L r ei, n re i) = (4, +1). Symbols show the numerically deter- 
mined energies while dashed lines show the energies predicted 
by the effective dimer-dimer model (see text for details). 



ural parity states are effected by the zero-range interac- 
tions. 

Next, we discuss the energy spectrum in the strongly- 
interacting regime and in the a^ aa ^ — > + limit. For 
small cS aa \ the energetically lowest lying states with 
even L Te \ belong to the "dimer plus dimer" family while 
the energetically lowest lying states with odd L rc \ be- 
long to the "dimer plus atom plus atom" family. Mo- 
tivated by this observation (see also Refs. [HI, [40|, HtJ) 
and by our discussion of the three-fermion system (see 
Sec. IIII A[) . Figs. ITUTa). (b) and (c) show the scaled en- 
ergies £W(2, 2) - 2i? d i mcr for L re i = 0,2 and 4. The 



rel 



FIG. 11: (Color online) Symbols show the difference between 
the lowest coefficient A" un it obtained from our numerically 
determined four-fermion energies and the lowest coefficient 
-Kunit, model predicted by the simple dimer-dimer model dis- 
cussed in the text as a function of L rc i (L re i even). To guide 
the eye symbols are connected by dotted lines. 



scaled four-fermion spectra for L IC \ — 2 and 4 con- 
tain a set of nearly constant scaled energies given by 
£ rcl (2,2) - 2£ dimcr fa (3/2 + L rcl + 2q)hw, where q = 
0,1,- ■■ [see dashed lines in Figs. ITUTb) and (c)]. Fig- 
ures [TUFb) and (c) show that this description improves 
with increasing L rel . The set of constant scaled ener- 
gies E re \(2, 2) — 2_Edimor is predicted to exist within the 
effective two-boson model, which treats the composite 
bosons as non-interacting if L rc \ > and L rel even [see 
the discussion below Eq. (fT5)l in Sec. [II]. Using that a 
subset of the scaled energies is approximately constant 
and that the dimer energy -Edimcr equals hu/2 at uni- 
tarity, the lowest -fT un it coefficient is, for L le \ even and 
L re i > 0, approximately given by X uni t, mo dci = L Tei + 1- 
The energy levels associated with these K un i t coefficients 
have the minimally allowed number of excitations in the 
hyperangular degrees of freedom and q excitations along 
the hyperradial coordinate R; these states transform, as 
a {aa) conges f rom oo to + , to "dimer plus dimer" states 
with 7i c ff radial excitations. Since the dimers are compos- 
ite bosons, an analogous set of states does not exist for 
odd L re i- Squares in Fig. [TT] show the difference between 
the lowest K un - lt coefficient and K un it, model, where i^ un it 
is taken from Tables fVll and IVlTl for L rc \ — 2 — 8 and L re i 
even. This figure supports our conjecture that the lowest 
coefficient K unlt converges to the value i^ U nit, model in the 
large L re \ limit. 

The behavior of the scaled energies for L xe \ = is 
different from that of the scaled energies for L IC \ > 
(L ro i even). Figure [TOf a) does not show a set of nearly 
constant scaled energy levels but instead shows a set of 
scaled energies that drop by approximately fiw in going 
from a( aa ) — > 0~ to cS aa ^ = oo and by another ftw in 



going from a 



(aa) 



OO to CL 



(aa) 



+ . Dashed lines 



in Fig. [TU] show the energies predicted by the effective 
dimer-dimer model. To obtain these energies, we solve 
Eq. (fTS]) for k = dd and use Eq. ([T4"]) to express the 
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effective dimer-dimer scattering length 

a (dd) in 

terms of 

the atom-atom scattering length a^ aa ^ throughout the en- 
tire crossover. The dashed lines agree surprisingly well 
with the full four-fermion energies throughout the en- 
tire crossover. Comparison with Fig. [Ha) shows that 
the effective two-particle model provides a better de- 
scription of the four-fermion than of the three-fermion 
system, particularly for negative a < -- aa K Intuitively, this 
might be explained by the fact that the three-fermion 
system contains an unpaired atom while all atoms par- 
ticipate in the molecule formation in the four-fermion 
system. Figure QJJ] suggests that the four-fermion en- 
ergy spectrum can be interpreted by considering avoided 
crossings between dimer-dimer states and states that be- 
long to the "dimer plus atom plus atom" and the "gas- 
like state" families. An analysis along these lines has 
been performed in Refs. [1(1113]; however, to the best of 
our knowledge this type of analysis has not been previ- 
ously applied to estimate the value of the lowest K un n 
coefficient for L re i = 0. Using that the energy in the 
aH 0+ limit is given by 2K dimor + (3/2 + 2n cS )haj 
and that the energy at unitarity is — using the argument 
above — 2E<n meI + (5/2 + 2q)fkv, the lowest -K un it coeffi- 
cient for L rc i = can be estimated to be -Kunit, model = 2. 
Table [VTI shows that this model is surprisingly accurate, 
i.e., K unit = 2.009 (see also Fig. ITT]) . 

We find numerically that the second lowest -Kunit co- 
efficient for even L Te \ (H Te ] = +1) and the lowest -K un it 
coefficient for odd L IC \ (II rc i = —1) appear to approach 
L le i + 7/4 for large L le \. At present, we have no expla- 
nation for this observation. While our analysis presented 
above allows for the prediction of the lowest -K un it co- 
efficient for all even L Te \, the remaining coefficients at 
unitarity appear to be result of an intricate interplay of 
how to best distribute the angular momentum among the 
six interparticle distances. 



IV. SUMMARY 

This paper considered s-wave interacting three- and 
four-fermion systems under external spherically symmet- 
ric harmonic confinement. If the range ro of the underly- 
ing interaction potential is much smaller than the other 
length scales of the system, i.e., the atom-atom s-wave 
scattering length a,( aa ) and the harmonic oscillator length 
a no , then the physics of two-component equal- mass Fermi 
systems is universal. If the s-wave scattering length is 
tuned in the vicinity of a so-called broad Feshbach reso- 
nance [EElfiSlzSli then the low-energy properties of these 
few-fermion systems are expected to behave universally, 
i.e., independently of the details of the underlying inter- 
action potential. Our results are expected to apply to 
the experimentally most frequently studied species 6 Li 
and 40 K [UllMi- 

We have determined the zero-range energy spectrum 
of the three- and four-fermion systems numerically for a 
wide range of s-wave scattering lengths a t - aa \ Our study 



of the three-fermion system is based on the Lippmann- 
Schwinger equation, which reduces the problem to solv- 
ing a set of coupled equations for each angular mo- 
mentum L re i [33|]. The three-fermion spectrum is ana- 
lyzed and interpreted following a number of different ap- 
proaches. In the weakly-interacting regime, the interac- 
tions are treated as a perturbation to the non-interacting 
trapped atomic Fermi gas. This approach correctly de- 
scribes the gas- like states of the system for small |a( aa )|, 
with cS aa ^ positive and negative, but does not describe 
states that are associated with the formation of diatomic 
molecules. The energies of this family of states are 
described by an effective two-particle Hamiltonian that 
treats the diatomic molecule as a composite point par- 
ticle and assumes that the spare atom and the molecule 
interact through the s-wave scattering channel character- 
ized by the effective atom-dimer scattering length a^ ad K 
Somewhat surprisingly, this effective two-particle model 
describes the energies of a subset of states not only quan- 
titatively correctly in the weakly-repulsive regime but 
also qualitatively for negative a^ aa ^ (see dashed lines in 
Fig. O. Building on this observation, we developed a 
simple model which predicts the -K un it coefficients with 
reasonably good accuracy. In doing so, our motivation 
was to develop a physical picture of the general features 
of the three-fermion energy spectrum that can, at least 
partially, be generalized to larger systems. A key re- 
sult of our analysis is that the energy levels determined 
by the lowest -K un it coefficient for each L rc i transform to 
"atom plus dimer" states as o/ aa ) changes from 00 to + . 
This "transformation" of the states can be interpreted 
nicely within the hyperspherical framework. Using hy- 
perspherical coordinates, the states for 

a (aa) = 

and 

+ are separable and the widths of the avoided crossings 
of the energy levels for finite are determined by the 
couplings between different hyperradial potential curves 
(see, e.g., Refs. (HlzzHzi as well as Ref. 

We also solved the time-independent Schrodinger equa- 
tion for the four-fermion system for a Gaussian two-body 
potential with varying range ro using the stochastic vari- 
ational approch. The resulting finite-range energies were 
then extrapolated to the ro — > limit. The energy spec- 
trum of the four-fermion system is much denser than 
that of the three-fermion system. Unlike for the three- 
fermion system, natural and unnatural parity states of 
the four-fermion system are shifted by the zero-range in- 
teractions. Our primary focus in this paper has been to 
characterize the energies of the lowest few states with 
natural parity throughout the crossover, including the 
infinitely strongly-interacting unitary regime. As in the 
three-fermion case, our semi- analytical perturbative ap- 
proach, which utilizes hyperspherical coordinates, repro- 
duces the numerically determined four-fermion energies 
with good accuracy in the weakly-attractive and weakly- 
repulsive regimes. We paid special attention to the 
infinitely strongly-interacting regime and analyzed the 
energy spectrum at unitarity within an effective two- 
boson model. This analysis provides a physical picture of 
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how the lowest i^" un it coefficient for natural parity states 
and even L rc \ comes about and how the associated lad- 
der of states transforms to "dimer plus dimer" states 
as a {aa ^ -> + . Furthermore, the fact that the effec- 
tive two-boson model reproduces a subset of energy lev- 
els semi-quantitatively throughout the entire crossover 
regime suggests that the concept of the effective dimer- 
dimer scattering length extends beyond the small 

a M regime (a< aa ) > 0). 

One motivation for studying small few-fermion sys- 
tems is to develop a "bottom-up approach" that investi- 
gates the microscopic physics of two-component Fermi 
gases by treating successively larger systems. With 
this motivation in mind, the three-fermion system can 



be considered the smallest system that models spin- 
imbalanced Fermi gases [33| while the four-fermion sys- 
tem can be considered the smallest system that models 
pairing physics throughout the BCS-BEC crossover of 
spin- balanced Fermi gases [40| . In the future, it will be 
interesting to extend the studies presented here to larger 
equal-mass few-fermion systems as well as to unequal- 
mass fermionic and bosonic systems. We also hope to 
use the perturbative treatment of the four-fermion sys- 
tem to estimate the fourth order virial coefficient in the 
weakly-interacting regimes. 

Support by the NSF through grant PHY-0855332 and 
by the ARO are gratefully acknowledged. 
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Auxiliary material for "Energy spectrum of 
harmonically trapped two-component Fermi gases: 
Three- and Four-Particle Problem" 

This EPAPS material contains the extrapolated zero- 
range energies E YC \(2,2) for negative scattering lengths 
a,( aa \ and the extrapolated scaled zero-range energies 
i? rc i(2,2) - 2i? dimcr and £' rol (2,2) - E dimeT for positive 
scattering lengths a*- 00 ) . The notation follows that of the 
main article. 
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TABLE VIII: Energies £ re i(2,2) for various L Icl and a (aa) . 
The energies (second through 8th column) are reported in 
units of Eho- Columns two through six show the energies 
for the five lowest states with (L re i,n rc i) = (0, +1) while 
columns seven through nine show the energies for the three 
lowest states with (L re i,II re i) = (1,-1). The uncertainty of 
the energies is estimated to be in the second to last digit for 
large |a* aa '| and 1 in the last digit reported for small |a' aa '|. 
"NF stands for non- interacting (a (aa) = 0). 
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L Te l = 1 
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5 


9594 


7.8830 


7.9276 


8 


1875 


8 


2514 


7 


0998 


7.1379 


7 


2671 


-4 


5 


8989 


7.8159 


7.8658 


8 


1528 


8 


2216 


7 


0571 


7.0983 


7 


2413 


-7/2 


5 


8238 


7.7337 


7.7899 


8 


1097 


8 


1837 


7 


0046 


7.0495 


7 


2092 


-3 


5 


7284 


7.6308 


7.6948 


8 


0553 


8 


1342 


6 


9389 


6.9879 


7 


1687 


-5/2 


5 


6040 


7.4994 


7.5732 


7 


9845 


8 


0667 


6 


8544 


6.9085 


7 


1158 


-2 


5 


4361 


7.3270 


7.4133 


7 


8899 


7 


9700 


6 


7432 


6.8027 


7 


0448 


-3/2 


5 


2009 


7.0936 


7.1990 


7 


7588 


7 


8230 


6 


5916 


6.6573 


6 


9456 


-1 


4 


8569 


6.7641 


6.9062 


7 


5698 


7 


5847 


6 


3779 


6.4497 


6 


8009 


-1/2 


4 


3324 


6.2745 


6.5041 


7 


1829 


7 


2863 


6 


0663 


6.1426 


6 


5796 







5092 


5.5101 


5.9441 


6 


5290 


6 


8482 


5 


5978 


5.6758 


6 


2305 
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TABLE IX: Nine lowest energies E iei (2, 2) for (L re i,n re i) = 
(2, +1) and various a ( - aa \ The energies (second through 10th 
column) are reported in units of E^ . The uncertainty of the 
energies is estimated to be in the second to last digit for large 
|a (aa) j and 1 in the last digit reported for small |a (aa) |. 



a ho /a 



(aa) 



NI 


13/2 


17/2 


17/2 


17/2 


17/2 


17/2 


17/2 


17/2 


17/2 


— 100 


6 


4801 


8.4766 


8.4789 


8.4816 


8.4872 


8.4872 


8.4924 


8.4942 


8 


5001 


—90 


6 


4779 


8.4740 


8.4765 


8.4796 


8.4858 


8.4858 


8.4915 


8.4936 


8 


5001 


—80 


6 


4751 


8.4706 


8.4735 


8.4769 


8.4838 


8.4838 


8.4903 


8.4927 


8 


5001 


—70 


6 


4716 


8.4665 


8.4698 


8.4737 


8.4817 


8.4817 


8.4891 


8.4917 


8 


5001 


—60 


6 


4669 


8.4610 


8.4648 


8.4694 


8.4787 


8.4787 


o a a fro 

8.4873 


8.4903 


8 


5001 


—50 


6 


4604 


8.4532 


8.4579 


8.4633 


8.4745 


8.4745 


8.4847 


8.4884 


8 


5001 


—40 


6 


4506 


8.4416 


8.4474 


8.4541 


8.4681 


8.4682 


8.4809 


8.4855 


8 


5002 


—30 


6 


4344 


8.4223 


8.4300 


8.4388 


8.4576 


8.4578 


8.4746 


8.4807 


8 


5002 


—20 


6 


4023 


8.3841 


8.3957 


8.4085 


8.4367 


8.4372 


8.4619 


8.4711 


8 


5003 


— 18 


6 


3918 


8.3714 


8.3843 


8.3984 


8.4298 


8.4304 


8.4577 


8.4679 


8 


5003 


— 16 


6 


3786 


8.3556 


8.3703 


8.3858 


8.4213 


8.4220 


8.4524 


8.4639 


8 


5004 


— 14 


6 


3619 


8.3354 


8.3522 


8.3697 


8.4104 


8.4112 


8.4457 


8.4588 


8 


5004 


— 12 


6 


3398 


8.3087 


8.3285 


a O A O A 

8.3484 


8.3959 


8.3970 


8.4367 


8.4520 


8 


5005 


— 10 


6 


3094 


8.2717 


8.2956 


8.3187 


8.3760 


8.3774 


8.4243 


8.4425 


8 


5005 


— 19/2 


6 


2992 


8.2579 


8.2841 


8.3079 


8.3695 


8.3712 


8.4205 


8.4394 


8 


5012 


—9 


6 


2887 


8.2450 


8.2728 


8.2976 


8.3626 


8.3645 


8.4161 


8.4361 


8 


5013 


— 17/2 


6 


2771 


8.2306 


8.2602 


8.2861 


8.3550 


8.3570 


8.4113 


8.4324 


8 


5013 


—8 


6 


2642 


8.2145 


8.2462 


8.2732 


8.3464 


8.3487 


8.4058 


8.4283 


8 


5014 


— 15/2 


6 


2494 


8.1964 


8.2304 


8.2587 


8.3368 


8.3393 


8.3997 


8.4235 


8 


5014 


—7 


6 


2329 


8.1759 


8.2127 


8.2422 


8.3260 


8.3287 


8.3927 


8.4182 


8 


5015 


— 13/2 


6 


2141 


8.1525 


8.1924 


8.2234 


8.3136 


8.3166 


8.3848 


8.4120 


8 


5016 


—6 


6 


1925 


8.1254 


8.1692 


8.2018 


8.2994 


8.3028 


8.3755 


8.4049 


8 


5016 


— 11/2 


6 


1675 


8.0939 


8.1423 


8.1765 


8.2828 


8.2867 


8.3646 


8.3965 


8 


5017 


-5 


6 


1380 


8.0566 


8.1109 


8.1467 


8.2634 


8.2678 


8.3518 


8.3864 


8 


5018 


-9/2 


6 


1030 


8.0121 


8.0735 


8.1112 


8.2402 


8.2453 


8.3362 


8.3743 


8 


5020 


-4 


6 


0607 


7.9581 


8.0287 


8.0680 


8.2122 


8.2181 


8.3172 


8.3593 


8 


5021 


-7/2 


6 


0087 


7.8911 


7.9738 


8.0147 


8.1777 


8.1847 


8.2933 


8.3402 


8 


5022 


-3 


5 


9432 


7.8064 


7.9053 


7.9475 


8.1343 


8.1428 


8.2627 


8.3155 


8 


5024 


-5/2 


5 


8587 


7.6964 


7.8178 


7.8606 


8.0782 


8.0888 


8.2219 


8.2819 


8 


5025 


-2 


5 


7459 


7.5496 


7.7027 


7.7451 


8.0035 


8.0170 


8.1658 


8.2342 


8 


5028 


-3/2 


5 


5895 


7.3473 


7.5456 


7.5865 


7.9001 


7.9180 


8.0847 


8.1620 


8 


5030 


-1 


5 


3613 


7.0597 


7.3211 


7.3615 


7.7501 


7.7757 


7.9610 


8.0430 


8 


5033 


-1/2 


5 


0077 


6.6422 


6.9799 


7.0320 


7.5194 


7.5610 


7.7613 


7.8277 


8 


5036 





4 


4185 


6.0390 


6.4201 


6.5391 


7.1288 


7.2218 


7.4247 


7.4271 


8 


0416 
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TABLE X: Energies E m i(2, 2) for various L rc i and a (aa) . The 
energies (second through 5th column) are reported in units of 
Eho- Column two shows the energy for the lowest state with 
(L rB i,n re i) = (3,-1) while columns three through five show 
the energies for the three lowest states with (L re i,II re i) = 
(4, +1). The uncertainty of the energies is estimated to be in 
the second to last digit for large |a' aa '| and 1 in the last digit 
reported for small la' 00 ']. 



a ho /a (aa) 




irel = 4 


NI 


15/2 


17/2 


17/2 


17/2 


-100 


7.4861 


8.4823 


8.4901 


8.4908 


-90 


7.4845 


8.4804 


8.4890 


8.4898 


-80 


7.4826 


8.4778 


8.4874 


8.4884 


-70 


7.4802 


8.4748 


8.4858 


8.4869 


-60 


7.4769 


8.4706 


8.4835 


8.4847 


-50 


7.4723 


8.4648 


8.4802 


8.4817 


-40 


7.4655 


8.4561 


8.4753 


8.4771 


-30 


7.4541 


8.4416 


8.4671 


8.4695 


-20 


7.4318 


8.4130 


8.4509 


8.4544 


-18 


7.4244 


8.4036 


8.4456 


8.4494 


-16 


7.4153 


8.3918 


8.4389 


8.4432 


-14 


7.4036 


8.3769 


8.4304 


8.4352 


-12 


7.3882 


8.3571 


8.4192 


8.4247 


-10 


7.3671 


8.3298 


8.4036 


8.4100 


-19/2 


7.3602 


8.3206 


8.3987 


8.4053 


-9 


7.3529 


8.3112 


8.3933 


8.4002 


-17/2 


7.3449 


8.3007 


8.3873 


8.3945 


-8 


7.3359 


8.2890 


8.3807 


8.3881 


-15/2 


7.3258 


8.2759 


8.3732 


8.3810 


-7 


7.3144 


8.2610 


8.3646 


8.3728 


-13/2 


7.3014 


8.2440 


8.3549 


8.3636 


-6 


7.2865 


8.2245 


8.3437 


8.3529 


11/0 




8.2018 


8.3307 


8.3403 


-5 


7.2490 


8.1751 


8.3153 


8.3255 


-9/2 


7.2250 


8.1432 


8.2969 


8.3077 


-4 


7.1961 


8.1046 


8.2745 


8.2860 


-7/2 


7.1606 


8.0569 


8.2469 


8.2588 


-3 


7.1162 


7.9965 


8.2119 


8.2242 


-5/2 


7.0593 


7.9179 


8.1662 


8.1786 


-2 


6.9841 


7.8121 


8.1046 


8.1163 


-3/2 


6.8812 


7.6630 


8.0182 


8.0275 


-1 


6.7346 


7.4409 


7.8909 


7.8944 


-1/2 


6.5164 


7.0871 


7.6837 


7.6926 





6.1756 


6.4855 


7.3384 


7.3677 
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TABLE XI: Scaled energies E ie \(2, 2) - 2_E d imor for vari- 
ous L re i and a ( - aa \ The scaled energies (second to 10th col- 
umn) are reported in units of -Eho- Columns two through 
four show the scaled energies for the three lowest states with 
(Lrci,n ro i) = (0,+l), columns five through seven show the 
scaled energies for the three lowest states with (Lrci,n r ci) = 
(2, +1) and columns eight through ten show the scaled en- 
ergies for the three lowest states with (Lr e i,n re i) = (4, +1). 
The uncertainty of the scaled energies is estimated to be in 
the last digit reported. 



a h o/a (aa) 


Lrel = 


L re l — 2 






Lrol = 4 




1/2 


2 


281 


4.357 


5 


190 


3 


457 


5 


239 


5 


539 


5 


495 


6.836 


6.890 


1 


2 


082 


4.206 


5 


718 


3 


486 


5 


464 


6 


073 


5 


499 


7.493 


7.719 


3/2 


1 


939 


4.072 


6 


138 


3 


497 


5 


486 


7 


264 


5 


500 


7.503 


9.033 


2 


1 


843 


3.967 


6 


045 


3 


499 


5 


496 


7 


500 


5 


500 


7.501 


9.508 


5/2 


1 


778 


3.888 


5 


963 


3 


501 


5 


499 


7 


504 


5 


500 


7.501 


9.506 


3 


1 


732 


3.831 


5 


899 


3 


500 


5 


500 


7 


506 


5 


500 


7.501 


9.504 


7/2 


1 


699 


3.787 


5 


848 


3 


501 


5 


500 


7 


503 


5 


500 


7.501 


9.504 


4 


1 


673 


3.752 


5 


807 


3 


502 


5 


500 


7 


502 


5 


500 


7.501 


9.501 


9/2 


1 


653 


3.724 


5 


773 


3 


501 


5 


500 


7 


502 


5 


500 


7.501 


9.501 


5 


1 


637 


3.701 


5 


745 


3 


500 


5 


500 


7 


502 


5 


500 


7.501 


9.501 



TABLE XII: Scaled energies _E re i(2,2) — -Edimor for various 
Lrd and d'°°'. The scaled energies (second through 5th col- 
umn) are reported in units of Eho- Columns two through 
four show the scaled energies for the three lowest states with 
(Lrd, II r ei) = (1, —1) and column five shows the scaled energy 
for the lowest state with {L ve \,W Ta \) = (3,-1). The uncer- 
tainty of the scaled energies is estimated to be in the last 
digit reported. 



a ho /a (aa) 


Lrel = 1 


Lrel = 3 


1/2 


4.927 


4.999 


5.703 


5.667 


1 


4.774 


4.834 


5.708 


5.697 


3/2 


4.654 


4.699 


5.745 


5.753 


2 


4.562 


4.593 


5.795 


5.813 


5/2 


4.490 


4.511 


5.841 


5.863 


3 


4.433 


4.446 


5.880 


5.900 


7/2 


4.387 


4.395 


5.909 


5.926 


4 


4.350 


4.352 


5.931 


5.945 


9/2 


4.318 


4.318 


5.947 


5.958 


5 


4.289 


4.290 


5.959 


5.968 



